#=======================================================================================#
#
#  R code to replicate analyses reported in the paper titled
#  "Human Rights Promotion and Democratic Allies."
#  R version is 4.2.3.
#
#  Author: Yasuki Kudo
#
#  Last modified: 2024-07-02
#=======================================================================================#

rm(list=ls())

sink("outputfile.txt")

# Packages
#=======================================================================================#
library(plm)              # ver. 2.6.3
library(tidyverse)        # ver. 2.0.0
library(marginaleffects)  # ver. 0.11.0
library(modelsummary)     # ver. 1.3.0
library(statar)           # ver. 0.7.4
library(gridExtra)        # ver. 2.3
library(grid)             # ver. 4.2.3
library(haven)            # ver. 2.5.4
library(ggridges)         # ver. 0.5.4
library(interflex)        # ver. 1.2.6
library(lmtest)           # ver. 0.9.40
library(tibble)           # ver. 3.2.1

# Analyses
#=======================================================================================#
#### import data ####
load("data/kudo2024isq_data.rdata")

#### settings of figures ####
line <- scale_linetype_manual(values=c("solid", "dashed"), labels=c("0", "1"))
color <- scale_color_manual(values=c("gray50", "black"), labels=c("1", "0"))
shape <- scale_shape_manual(values=c(16, 15), labels=c("5th pctl", "95th pctl"))

#### main models ####
for(h in 1:5) {
  
  # h=1: main models (Table 1 and Figure 1)
  # h=2: models with all states (Appendix A5)
  # h=3: models with fewer variables (Appendix A3)
  # h=4: 40th pctl threshold (Appendix A2)
  # h=5: 60th pctl threshold (Appendix A2)
  
  if(h %in% c(1, 2, 3)) {
    df$dmpatron <- df$dmpatron5
    df$ndmpatron <- df$ndmpatron5
    df$dmally <- df$dmally5
    df$ndmally <- df$ndmally5
  }
  if(h == 4) {
    df$dmpatron <- df$dmpatron4
    df$ndmpatron <- df$ndmpatron4
    df$dmally <- df$dmally4
    df$ndmally <- df$ndmally4
  }
  if(h == 5) {
    df$dmpatron <- df$dmpatron6
    df$ndmpatron <- df$ndmpatron6
    df$dmally <- df$dmally6
    df$ndmally <- df$ndmally6
  }
  
  for(i in 1:4) {
    
    # i=1: Model 1
    # i=2: Model 2
    # i=3: Model 3
    # i=4: Model 4
    
    if(i %in% c(1, 2)) {df$hr <- df$MVAphyFariss}
    if(i %in% c(3, 4)) {df$hr <- df$MVAempFariss}
    
    if(h != 3) {
      
      if(i %in% c(1, 3)) {
        
        form <- hr ~ 
          dmpatron + yb + ndmpatron + dmally + ndmally +
          loggdppc + logpop + logoil + threats + logodagdp + logtradegdp + logarmsgdp + logtroops +
          empFariss + phyFariss + HRregime + multi + coldwar
        
      }
      
      if(i %in% c(2, 4)) {
        
        form <- hr ~ 
          dmpatron*yb + ndmpatron + dmally + ndmally +
          loggdppc + logpop + logoil + threats + logodagdp + logtradegdp + logarmsgdp + logtroops +
          empFariss + phyFariss + HRregime + multi + coldwar
        
      }
      
    }else{
      
      if(i %in% c(1, 3)) {
        
        form <- hr ~ 
          dmpatron + yb + loggdppc + logpop + threats + logoil + empFariss + phyFariss + multi + coldwar
        
      }
      
      if(i %in% c(2, 4)) {
        
        form <- hr ~ 
          dmpatron*yb + loggdppc + logpop + threats + logoil + empFariss + phyFariss + multi + coldwar
        
      }
    }
    
    if(h != 2) {
      
      reg <- plm(form, data = df, subset = autocracy == 1, effect = "individual", model = "pooling", index = c("cowcode","year"))
      
    }else{
      
      reg <- plm(form, data = df, effect = "individual", model = "pooling", index = c("cowcode","year"))
      
    }
    
    vcov.bk <- vcovBK(reg, cluster = "time", type = "HC0")
    mod <- coeftest(reg, vcov.bk, save = T)
    
    if(i %in% c(2, 4)) {
      
      if(h == 1 && i == 2) {title <- ggtitle("Model 2 (DV: Physical Integrity Rights)")}
      if(h == 1 && i == 4) {title <- ggtitle("Model 4 (DV: Empowerment Rights)")}
      if(h == 2 && i == 2) {title <- ggtitle("Model A14 (DV: Physical Integrity Rights)")}
      if(h == 2 && i == 4) {title <- ggtitle("Model A16 (DV: Empowerment Rights)")}
      if(h == 3 && i == 2) {title <- ggtitle("Model A10 (DV: Physical Integrity Rights)")}
      if(h == 3 && i == 4) {title <- ggtitle("Model A12 (DV: Empowerment Rights)")}
      if(h == 4 && i == 2) {title <- ggtitle("Model A2 (DV: Physical Integrity Rights)")}
      if(h == 4 && i == 4) {title <- ggtitle("Model A4 (DV: Empowerment Rights)")}
      if(h == 5 && i == 2) {title <- ggtitle("Model A6 (DV: Physical Integrity Rights)")}
      if(h == 5 && i == 4) {title <- ggtitle("Model A8 (DV: Empowerment Rights)")}
      if(h == 5 && i == 2) {title <- ggtitle("Model A6 (DV: Physical Integrity Rights)")}
      if(h == 5 && i == 4) {title <- ggtitle("Model A8 (DV: Empowerment Rights)")}
      
      if(h != 2) {yb.minmax <- quantile(subset(df, autocracy == 1)$yb, prob = seq(0, 1, length = 15), na.rm = T)}
      if(h == 2) {yb.minmax <- quantile(df$yb, prob = seq(0, 1, length = 15), na.rm = T)}
      
      comp <- comparisons(model = reg, vcov = vcov.bk, variables = list(dmpatron = c(0, 1)), 
                          newdata = datagrid(yb = yb.minmax), conf_level = 0.95)
      
      comp %>% 
        ggplot() +
        theme_bw() +
        geom_ribbon(aes(x = yb, ymax = conf.high, ymin = conf.low), alpha = 0.3) +
        geom_line(aes(x = yb, y = estimate), size = 0.5) +
        geom_hline(yintercept = 0, linetype = "dotted", size = 0.25) +
        geom_vline(xintercept = pctile(subset(df, autocracy == 1)$yb, na.rm = T, probs = c(0.05, 0.5, 0.95)), linetype = "dotted", size = 0.25) +
        theme(text = element_text(size = 9),
              plot.title = element_text(size = 9),
              panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank(),
              axis.text.x = element_text(size = 6,
                                         color = "black"),
              axis.text.y = element_text(size = 6,
                                         color = "black"),
              plot.margin=unit(c(3,3,3,3), 'mm')) +
        xlab("Youth Bulges") + ylab("Marginal Effect (95% CI)") + title ->> margplot
      
      if(h == 1) {
        
        pred <- predictions(model = reg, vcov = vcov.bk, variables = list(yb = yb.minmax), conf_level = 0.95,
                            newdata = datagrid(dmpatron = c(0, 1), ndmpatron = 0, dmally = 0, ndmally = 1,
                                               multi = 0, coldwar = 0))
        
        if(i == 2) {
          theme <- 
            theme(text = element_text(size = 9),
                  plot.title = element_text(size = 9),
                  strip.text = element_text(size = 9),
                  strip.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
                  panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                  axis.text.x = element_text(size = 6, color = "black"),
                  axis.text.y = element_text(size = 6, color = "black"),
                  legend.position=c(.74,.8),
                  legend.direction="vertical",
                  legend.key.width = unit(.6,"cm"),
                  legend.key.height = unit(.2,"cm"),
                  legend.text = element_text(size = 8),
                  legend.title = element_text(size = 8),
                  legend.box.background = element_rect(colour = "black"),
                  plot.margin=unit(c(3,3,3,3), 'mm'))
          
        }
        
        if(i == 4) {
          
          theme <- 
            theme(text = element_text(size = 9),
                  plot.title = element_text(size = 9),
                  strip.text = element_text(size = 9),
                  strip.background = element_rect(colour = 'black', fill = 'white', 
                                                  linetype='solid'),
                  panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                  axis.text.x = element_text(size = 6, color = "black"),
                  axis.text.y = element_text(size = 6, color = "black"),
                  legend.position="none",
                  plot.margin=unit(c(3,3,3,3), 'mm')) 
          
        }
        
        pred %>% 
          ggplot() +
          theme_bw() +
          line +
          geom_line(aes(x=yb, y=estimate, linetype = as.factor(dmpatron)), size = 0.5) +
          geom_vline(xintercept = pctile(subset(df, autocracy == 1)$yb, na.rm = T, probs = c(0.05, 0.5, 0.95)),
                     linetype = "dotted", size = 0.25) +
          theme +
          guides(linetype = guide_legend(ncol = 1, label.position = "right")) +
          xlab("Youth Bulges") + ylab("Predicted Value") +
          title +
          labs(linetype = "Democratic Patron") -> predplot
        
      }
    }
    
    if(h == 1 & i == 1) {
      mod ->> mod1.th5
      summary(reg)$r.squared[1] ->> r2.mod1.th5
    }
    
    if(h == 1 & i == 2) {
      mod ->> mod2.th5
      margplot ->> margplot2.th5
      predplot ->> predplot2.th5
      summary(reg)$r.squared[1] ->> r2.mod2.th5
    }
    
    if(h == 1 & i == 3) {
      mod ->> mod3.th5
      summary(reg)$r.squared[1] ->> r2.mod3.th5
    }
    
    if(h == 1 & i == 4) {
      mod ->> mod4.th5
      margplot ->> margplot4.th5
      predplot ->> predplot4.th5
      summary(reg)$r.squared[1] ->> r2.mod4.th5
    }
    
    if(h == 2 & i == 1) {
      mod ->> mod1.all
      summary(reg)$r.squared[1] ->> r2.mod1.all
    }
    
    if(h == 2 & i == 2) {
      mod ->> mod2.all
      margplot ->> margplot2.all
      summary(reg)$r.squared[1] ->> r2.mod2.all
    }
    
    if(h == 2 & i == 3) {
      mod ->> mod3.all
      summary(reg)$r.squared[1] ->> r2.mod3.all
    }
    
    if(h == 2 & i == 4) {
      mod ->> mod4.all
      margplot ->> margplot4.all
      summary(reg)$r.squared[1] ->> r2.mod4.all
    }
    
    if(h == 3 & i == 1) {
      mod ->> mod1.nested
      summary(reg)$r.squared[1] ->> r2.mod1.nested
    }
    
    if(h == 3 & i == 2) {
      mod ->> mod2.nested
      margplot ->> margplot2.nested
      summary(reg)$r.squared[1] ->> r2.mod2.nested
    }
    
    if(h == 3 & i == 3) {
      mod ->> mod3.nested
      summary(reg)$r.squared[1] ->> r2.mod3.nested
    }
    
    if(h == 3 & i == 4) {
      mod ->> mod4.nested
      margplot ->> margplot4.nested
      summary(reg)$r.squared[1] ->> r2.mod4.nested
    }
    
    if(h == 4 & i == 1) {
      mod ->> mod1.th4
      summary(reg)$r.squared[1] ->> r2.mod1.th4
    }
    
    if(h == 4 & i == 2) {
      mod ->> mod2.th4
      margplot ->> margplot2.th4
      summary(reg)$r.squared[1] ->> r2.mod2.th4
    }
    
    if(h == 4 & i == 3) {
      mod ->> mod3.th4
      summary(reg)$r.squared[1] ->> r2.mod3.th4
    }
    
    if(h == 4 & i == 4) {
      mod ->> mod4.th4
      margplot ->> margplot4.th4
      summary(reg)$r.squared[1] ->> r2.mod4.th4
    }
    
    if(h == 5 & i == 1) {
      mod ->> mod1.th6
      summary(reg)$r.squared[1] ->> r2.mod1.th6
    }
    
    if(h == 5 & i == 2) {
      mod ->> mod2.th6
      margplot ->> margplot2.th6
      summary(reg)$r.squared[1] ->> r2.mod2.th6
    }
    
    if(h == 5 & i == 3) {
      mod ->> mod3.th6
      summary(reg)$r.squared[1] ->> r2.mod3.th6
    }
    
    if(h == 5 & i == 4) {
      mod ->> mod4.th6
      margplot ->> margplot4.th6
      summary(reg)$r.squared[1] ->> r2.mod4.th6
    }
  }
}

#### long-term effects ####
comp.lt <- tibble()

for(i in 1:10){
  
  if(i == 1) {df$hr <- df$LD1phyFariss}
  if(i == 2) {df$hr <- df$LD2phyFariss}
  if(i == 3) {df$hr <- df$LD3phyFariss}
  if(i == 4) {df$hr <- df$LD4phyFariss}
  if(i == 5) {df$hr <- df$LD5phyFariss}
  if(i == 6) {df$hr <- df$LD1empFariss}
  if(i == 7) {df$hr <- df$LD2empFariss}
  if(i == 8) {df$hr <- df$LD3empFariss}
  if(i == 9) {df$hr <- df$LD4empFariss}
  if(i == 10) {df$hr <- df$LD5empFariss}
  
  form <- hr ~ 
    dmpatron5*yb + ndmpatron5 + dmally5 + ndmally5 +
    loggdppc + logpop + logoil + threats + logodagdp + logtradegdp + logarmsgdp + logtroops +
    empFariss + phyFariss + HRregime + multi + coldwar
  
  reg <- plm(form, data = df, subset = autocracy == 1, effect = "individual", model = "pooling", index = c("cowcode","year"))
  vcov.bk <- vcovBK(reg, cluster = "time", type = "HC0")
  new.data = datagrid(model = reg, yb = pctile(subset(df, autocracy == 1)$yb, na.rm = T, probs = c(0.05, 0.95)))
  comp <- comparisons(model = reg, vcov = vcov.bk, variables = list(dmpatron5 = c(0, 1)), newdata = new.data, conf_level = 0.95)
  comp.lt <- rbind(comp.lt, comp)
  
}

comp.lt.phy <- comp.lt[1:10,] %>%
  mutate(sig = ifelse(conf.low*conf.high > 0, "1", "0")) %>%
  mutate(t = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5))

comp.lt.emp <- comp.lt[11:20,] %>%
  mutate(sig = ifelse(conf.low*conf.high > 0, "1", "0")) %>%
  mutate(t = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5))

comp.lt.phy %>%
  ggplot() +
  theme_bw() +
  geom_pointrange(aes(x=as.factor(t), y=estimate, ymin = conf.low, ymax = conf.high,
                      shape = as.factor(yb), color = sig),
                  position = position_dodge(0.4),
                  size = .7, linewidth = .5) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.25) +
  color +　shape +
  theme(text = element_text(size = 9),
        plot.title = element_text(size = 9),
        strip.text = element_text(size = 9),
        strip.background = element_rect(colour = 'black',
                                        fill = 'white', 
                                        linetype='solid'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 8,
                                   color = "black"),
        axis.text.y = element_text(size = 6,
                                   color = "black"),
        legend.position="none",
        plot.margin=unit(c(3,3,3,3), 'mm')) +
  scale_x_discrete(labels= c("1" = "t+1", "2" = "t+2", "3" = "t+3", "4" = "t+4", "5" = "t+5")) +
  labs(shape = "Youth Bulges") +
  ylim(-0.27, 1) +
  xlab("") + ylab("Marginal Effect (95% CI)") +
  ggtitle("DV: Physical Integrity Rights") -> longplot.phy

comp.lt.emp %>%
  ggplot() +
  theme_bw() +
  geom_pointrange(aes(x=as.factor(t), y=estimate, ymin = conf.low, ymax = conf.high,
                      shape = as.factor(yb), color = sig),
                  position = position_dodge(0.4),
                  size = .7, linewidth = .5) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.25) +
  color + shape +
  theme(text = element_text(size = 9),
        plot.title = element_text(size = 9),
        strip.text = element_text(size = 9),
        strip.background = element_rect(colour = 'black',
                                        fill = 'white', 
                                        linetype='solid'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 8, color = "black"),
        axis.text.y = element_text(size = 6, color = "black"),
        legend.position=c(.77,.82),
        legend.direction="vertical",
        legend.box = "horizontal",
        legend.key.width = unit(.8,"cm"),
        legend.key.height = unit(.28,"cm"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.box.background = element_rect(colour = "black"),
        plot.margin=unit(c(3,3,3,3), 'mm')) +
  scale_x_discrete(labels= c("1" = "t+1", "2" = "t+2", "3" = "t+3", "4" = "t+4", "5" = "t+5")) +
  guides(color = F) +
  labs(shape = "Youth Bulges") +
  ylim(-0.27, 1) +
  xlab("") + ylab("Marginal Effect (95% CI)")  + 
  ggtitle("DV: Empowerment Rights") -> longplot.emp


#### binning estimators ####
yb_dmpt1 <-  subset(df, autocracy == 1 & dmpatron5 == 1)$yb
yb_dmpt0 <-  subset(df, autocracy == 1 & dmpatron5 == 0)$yb
yb_dens_dmpt1 <- approxfun(density(na.omit(yb_dmpt1)))
yb_dens_dmpt0 <- approxfun(density(na.omit(yb_dmpt0)))
df_auto <- subset(df, autocracy == 1) %>% 
  mutate(yb_dens = ifelse(dmpatron5 == 1, yb_dens_dmpt1(yb), yb_dens_dmpt0(yb)))

for(i in 1:2) {
  
  df$hr <- NA
  
  if(i == 1) {
    
    df$hr <- df$MVAphyFariss
    
  } else {
    
    df$hr <- df$MVAempFariss
    
  }
  
  y <- "hr"
  x <- "yb"
  d <- "dmpatron5"
  z <- c("ndmpatron5", "dmally5", "ndmally5", "loggdppc", "logpop",
         "logoil", "threats", "logodagdp", "logtradegdp", "logarmsgdp", "logtroops",
         "empFariss", "phyFariss", "HRregime", "multi", "coldwar")
  
  intf <- 
    interflex(Y = y, D = d, X = x, Z = z, estimator = "binning", na.rm = T, 
              vcov.type = "robust", nbins = 3, data = as.data.frame(subset(df, autocracy == 1)))
  
  if(i == 1) {
    
    intf[["est.bin"]][["1"]] ->> bin2
    intf[["est.lin"]][["1"]] ->> marg2
    
  }
  
  if(i == 2) {
    
    intf[["est.bin"]][["1"]] ->> bin4
    intf[["est.lin"]][["1"]] ->> marg4
    
  }
  
  if(i == 1) {
    
    title <- ggtitle("DV: Physical Integrity Rights \n (based on Model 2)")
    
  }
  
  if(i == 2) {
    
    title <- ggtitle("DV: Empowerment Rights \n (based on Model 4)")
    
  }
  
  if(i == 1) {
    
    marg <- marg2
    bin <- bin2
    
  } else {
    
    marg <- marg4
    bin <- bin4
    
  }
  
  ggplot(NULL) +
    theme_bw() +
    geom_hline(yintercept = 0, linetype = "dotted", size = .25) +
    geom_vline(xintercept = pctile(subset(df, autocracy == 1)$yb, na.rm = T,
                                   probs = c(0.05, 0.5, 0.95)),
               linetype = "dotted", size = .25) +
    geom_ribbon(aes(x=X, ymax =`upper CI(95%)`, ymin = `lower CI(95%)`),
                alpha = 0.3, data = as.data.frame(marg)) +
    geom_line(aes(x=X, y=TE), size = 0.5, data = as.data.frame(marg)) +
    geom_pointrange(aes(x = x0, y = coef, ymax = `CI.lower`, ymin = `CI.upper`),
                    size = 0.5, linewidth = 0.5, color = "red", data = as.data.frame(bin)) +
    geom_ridgeline(aes(x = yb, y = -0.4,  height=yb_dens/50, color = as.factor(dmpatron5)),
                   alpha = 0.3, size = 0.3, data = df_auto) +
    theme(text = element_text(size = 9),
          legend.position = "none",
          plot.title = element_text(size = 9),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          axis.text.x = element_text(size = 6,color = "black"),
          axis.text.y = element_text(size = 6, color = "black"),
          plot.margin=unit(c(3,3,3,3), 'mm')) +
    xlab("Youth Bulges") + ylab("Marginal Effect (95% CI)") + title -> binplot
  
  if(i == 1) {
    
    binplot ->> binplot2
    
  } else {
    
    binplot ->> binplot4
    
  }
}

#### additional analysis using UPR recommendations ####
TB2021 %>%
  ggplot(aes(yb, hr)) +
  theme_bw() +
  geom_point(aes(shape = type), size = 1.9, color = "gray50") +
  geom_smooth(aes(linetype = type), method = "lm", se = F, color = "black", size = .7) +
  scale_colour_identity(guide = "legend") +
  theme(text = element_text(size = 9),
        plot.title = element_text(size = 9),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 6, color = "black"),
        axis.text.y = element_text(size = 6, color = "black"),
        plot.margin=unit(c(3,7,3,3), 'mm'),
        legend.position=c(.67,.93),
        legend.direction="horizontal",
        legend.key.width = unit(0.6,"cm"),
        legend.key.height = unit(.15,"cm"),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.box.background = element_rect(colour = "black")) +
  scale_shape_discrete(name="") +
  scale_linetype_discrete(name="") +
  xlab("Youth Bulges") + ylab("Human Rights Score") -> hr_yb

TB2021 %>%
  filter(type == "PIR") %>%
  ggplot(aes(yb, mean_severity)) + 
  theme_bw() +
  geom_point(size = 1.9, color = "gray50") +
  geom_smooth(method = lm, se = F, color = "black", size = .7) +
  scale_colour_identity(guide = "legend") +
  theme(text = element_text(size = 9),
        plot.title = element_text(size = 9),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(size = 6, color = "black"),
        axis.text.y = element_text(size = 6, color = "black"),
        plot.margin=unit(c(3,7,3,3), 'mm')) +
  xlab("Youth Bulges") + ylab("Mean Severity of Recommendations")-> severity_yb

# Correlation
cor.test(subset(TB2021, type == "PIR")$yb, subset(TB2021, type == "PIR")$mean_severity, use = "complete.obs")
cor.test(subset(TB2021, type == "PIR")$yb, subset(TB2021, type == "PIR")$hr, use = "complete.obs")
cor.test(subset(TB2021, type == "ER")$yb, subset(TB2021, type == "ER")$hr, use = "complete.obs")


#### histograms ####
for(i in 1:7) {
  
  if(i == 1) {
    df_auto$hist.var <- df_auto$empFariss
    xlab <- xlab("Empowerment rights")
  }
  
  if(i == 2) {
    df_auto$hist.var <- df_auto$phyFariss
    xlab <- xlab("Physical integrity rights")
  }
  
  if(i == 3) {
    df_auto$hist.var <- df_auto$yb
    xlab <- xlab("Youth bulges")
  }
  
  if(i == 4) {
    df_auto$hist.var <- df_auto$dmpatron5
    xlab <- xlab("Democratic patron")
  }
  
  if(i == 5) {
    df_auto$hist.var <- df_auto$dmally5
    xlab <- xlab("Democratic ally")
  }
  
  if(i == 6) {
    df_auto$hist.var <- df_auto$ndmpatron5
    xlab <- xlab("Nondemocratic patron")
  }
  
  if(i == 7) {
    df_auto$hist.var <- df_auto$ndmally5
    xlab <- xlab("Nondemocratic ally")
  }
  
  if(i <= 3) {
    
    ggplot() +
      theme_bw() +
      geom_histogram(aes(hist.var), df_auto, color="black", fill="gray", bins=20, linewidth = .2) +
      theme(text = element_text(size = 8),
            strip.text = element_text(size = 8),
            strip.background = element_rect(colour = 'black',
                                            fill = 'white', 
                                            linetype='solid'),
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            axis.text.x = element_text(size = 6,
                                       color = "black"),
            axis.text.y = element_text(size = 6,
                                       color = "black")) +
      xlab + ylab("Frequency") -> hist.plot
    
  } else {
    
    ggplot() +
      theme_bw() +
      geom_bar(aes(as.factor(hist.var)), df_auto, color="black", fill="gray",
               stat="count", width=.5, linewidth = .2) +
      theme(text = element_text(size = 8),
            strip.text = element_text(size = 8),
            strip.background = element_rect(colour = 'black',
                                            fill = 'white', 
                                            linetype='solid'),
            panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            axis.text.x = element_text(size = 8,
                                       color = "black"),
            axis.text.y = element_text(size = 6,
                                       color = "black")) +
      xlab + ylab("Frequency") -> hist.plot
    
  }
  
  if(i == 1) hist.plot ->> histemp
  if(i == 2) hist.plot ->> histphy
  if(i == 3) hist.plot ->> histyb
  if(i == 4) hist.plot ->> histdmpatron
  if(i == 5) hist.plot ->> histdmally
  if(i == 6) hist.plot ->> histndmpatron
  if(i == 7) hist.plot ->> histndmally
  
}

# Display tables
#=======================================================================================#
coef_map <- c(dmpatron = "Democratic patron",
              ndmpatron = "Nondemocratic patron",
              dmally = "Democratic ally",
              ndmally = "Nondemocratic ally",
              yb = "Youth bulges",
              `dmpatron:yb` = "Democratic patron. x YB",
              loggdppc = "GDP per capita",
              logpop = "Total population",
              logoil = "Oil rents",
              threats = "External threats",
              logodagdp = "ODA from DAC",
              logtradegdp = "Trade with democracies",
              logarmsgdp = "Arms imports from democracies",
              logtroops = "US troops",
              HRregime = "Human rights regime",
              multi = "NATO, OAS, or AU",
              coldwar = "Cold War",
              empFariss = "Empowerment rights",
              phyFariss = "Physical integrity rights")

# Table 1
rows.t1 <- tribble(~term, ~`Model 1`,  ~`Model 2`, ~`Model 3`, ~`Model 4`,
                   "R2", r2.mod1.th5, r2.mod2.th5, r2.mod3.th5, r2.mod4.th5)
modelsummary(list("Model 1" = mod1.th5, "Model 2" = mod2.th5,
                  "Model 3" = mod3.th5, "Model 4" = mod4.th5),
             align = "ldddd",
             stars = c('+' = .1, '*' = .05, '**' = .01), 
             coef_map = coef_map,
             gof_map = c("nobs"),
             title = "Table 1",
             add_rows = rows.t1)

# Table A1
rows.ta1 <- tribble(~term, ~`Model A1`,  ~`Model A2`, ~`Model A3`, ~`Model A4`,
                    "R2", r2.mod1.th4, r2.mod2.th4, r2.mod3.th4, r2.mod4.th4)
modelsummary(list("Model A1" = mod1.th4, "Model A2" = mod2.th4,
                  "Model A3" = mod3.th4, "Model A4" = mod4.th4),
             align = "ldddd",
             stars = c('+' = .1, '*' = .05, '**' = .01), 
             coef_map = coef_map,
             gof_map = c("nobs"),
             title = "Table A1 (40th pctl)",
             add_rows = rows.ta1)

# Table A2
rows.ta2 <- tribble(~term, ~`Model A5`,  ~`Model A6`, ~`Model A7`, ~`Model A8`,
                    "R2", r2.mod1.th6, r2.mod2.th6, r2.mod3.th6, r2.mod4.th6)
modelsummary(list("Model A5" = mod1.th6, "Model A6" = mod2.th6, 
                  "Model A7" = mod3.th6, "Model A8" = mod4.th6),
             align = "ldddd",
             stars = c('+' = .1, '*' = .05, '**' = .01), 
             coef_map = coef_map,
             gof_map = c("nobs"),
             title = "Table A2 (60th pctl)",
             add_rows = rows.ta2)

# Table A3
rows.ta3 <- tribble(~term, ~`Model 1`,  ~`Model 2`, ~`Model 3`, ~`Model 4`,
                    "R2", r2.mod1.nested, r2.mod2.nested, r2.mod3.nested, r2.mod4.nested)
modelsummary(list("Model A9" = mod1.nested, "Model A10" = mod2.nested,
                  "Model A11" = mod3.nested, "Model A12" = mod4.nested),
             align = "ldddd",
             stars = c('+' = .1, '*' = .05, '**' = .01), 
             coef_map = coef_map,
             gof_map = c("nobs"),
             title = "Table A3 (Fewer variables)",
             add_rows = rows.ta3)

# Table A4
rows.ta4 <- tribble(~term, ~`Model 1`,  ~`Model 2`, ~`Model 3`, ~`Model 4`,
                    "R2", r2.mod1.all, r2.mod2.all, r2.mod3.all, r2.mod4.all)
modelsummary(list("Model A13" = mod1.all, "Model A14" = mod2.all,
                  "Model A15" = mod3.all, "Model A16" = mod4.all),
             align = "ldddd",
             stars = c('+' = .1, '*' = .05, '**' = .01), 
             coef_map = coef_map,
             gof_map = c("nobs"),
             title = "Table A4 (All states)",
             add_rows = rows.ta4)

# Export figures
#=======================================================================================#

# Figure 1
grid.arrange(arrangeGrob(margplot2.th5, margplot4.th5,
                         top=textGrob("Effects of Democratic Patron", gp=gpar(fontsize=9)),
                         ncol=1),
             arrangeGrob(predplot2.th5, predplot4.th5,
                         top=textGrob("Predicted Human Rights Scores", gp=gpar(fontsize=9)),
                         ncol=1), ncol = 2) -> plotAuto
ggsave("figure1.jpg", plotAuto, width = 6, height = 5.5, units = "in", dpi = 2000)  

# Figure 2
grid.arrange(longplot.phy, longplot.emp, ncol = 2) -> lag_plots
ggsave("figure2.jpg", lag_plots, width = 6, height = 3, units = "in", dpi = 2500) 

# Figure A1
grid.arrange(histemp, histphy, histyb, histdmpatron, 
             histndmpatron, histdmally, histndmally, ncol = 4) -> hist
ggsave("figureA1.jpg", hist, width = 6, height = 4, units = "in", dpi = 2500) 

# Figure A2
grid.arrange(margplot2.th4, margplot4.th4, ncol = 2) -> plotAutoR4
ggsave("figureA2.jpg", plotAutoR4, width = 6, height = 3, units = "in", dpi = 2000)  

# Figure A3
grid.arrange(margplot2.th6, margplot4.th6, ncol = 2) -> plotAutoR6
ggsave("figureA3.jpg", plotAutoR6, width = 6, height = 3, units = "in", dpi = 2000)  

# Figure A4
grid.arrange(margplot2.nested, margplot4.nested, ncol = 2) -> plotNested
ggsave("figureA4.jpg", plotNested, width = 6, height = 3, units = "in", dpi = 2000)  

# Figure A5
grid.arrange(binplot2, binplot4, ncol = 2) -> plotBin
ggsave("figureA5.jpg", plotBin, width = 6, height = 3, units = "in", dpi = 2000)  

# Figure A6
grid.arrange(margplot2.all, margplot4.all, ncol = 2) -> plotAll
ggsave("figureA6.jpg", plotAll, width = 6, height = 3, units = "in", dpi = 2000)  

# Figure A7
grid.arrange(arrangeGrob(hr_yb, bottom=textGrob("Human Rights Score v. Youth Bulges", gp = gpar(fontsize=9)),
                         ncol=1),
             arrangeGrob(severity_yb, bottom=textGrob("Severity of Recomm. v. Youth Bulges", gp = gpar(fontsize=9)),
                         ncol=1),
             ncol = 2) -> cor_plot
ggsave("figureA7.jpg", cor_plot, width = 6, height = 3, units = "in", dpi = 2500)  


sink()

# [End of script]
